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Abstract 

A model of static boson-fermion star with spherical symmetry based on the scalar-tensor theory 
of gravity with massive dilaton field is investigated numerically. 

Since the radius of star is a priori an unknown quantity, the corresponding boundary value 
problem (BVP) is treated as a nonlinear spectral problem with a free internal boundary. The 
Continuous Analogue of Newton Method (CANM) for solving this problem is applied. 

Information about basic geometric functions and the functions describing the matter fields, 
which build the star is obtained. In a physical point of view the main result is that the structure 
and properties of the star in presence of massive dilaton field depend essentially both of its fermionic 
and bosonic components. 

Keywords, boson-fermion star, scalar-tensor theory of gravity, massive dilaton field, two- 
parametric nonlinear spectral problem, continuous analog of Newton method, method of spline- 
collocation. 

Subject classification: 65C20, 65P30, 83-08, 83D05. 



1 Introduction 

The most natural and promising generalizations of general relativity are the scalar-tensor theories of 
gravity ~ 0. In these theories gravity is mediated not only by a tensor field (the metric of space- 
time) but also by a scalar field (the dilaton). The scalar-tensor theories of gravity contain arbitrary 
functions of the scalar field that determine the gravitational "constant" as a dynamical variable and 
the strength of the coupling between the scalar field and matter. It should be stressed that specific 
scalar-tensor theories of gravity arise naturally as a low energy limit of the string theory ||] - ]l3t , 
which is the most promising modern model of the unification of all fundamental physical interactions. 

*This research was supported by the Bulgarian Ministry of Education, Science and Technologies under the Grants 
NoNo MM-602/96, F610/99 and by the Sofia University Research Fund, Contr. No 245/99. 
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If the string theory and its low energy limit are relevant to the real world, then the dilaton must be 
massive [p^ . Unfortunately, our current understanding of how the dilaton acquires mass is primitive 
and it is tied to our lack of understanding of supersymmetry breaking. At present, we do not have a 
model of how the dilaton mass is generated in the string theory. Besides the mass term for dilaton field 
we may consider the general case of arbitrary dilaton potential, describing its nonlinear self-interaction. 

From physical point of view, it is important to know how the dilaton mass and, in general, the 
dilaton potential influence the structure and stability of compact objects such as neutron stars, boson 
stars, and mixed fermion-boson stars. 

It is known that the predictions of scalar-tensor theories of gravity with massless dilaton may 
differ drastically from these of general relativity. For example, the phenomenon of "spontaneous 
scalarization" was discovered recently [|l5|, as a non-perturbative strong field effect in a massive 
neutron star. The existence of this effect poses some important physical questions That is why it 
is natural to ask whether or not the "spontaneous scalarization" will occur when the dilaton is massive. 
In recent years, the boson stars in scalar-tensor theories of gravity with massless dilaton have been 
widely studied both analytically and numerically (see for example jl^ - [^). The study of boson stars 
in the case of massive dilaton is physically interesting and may be important for the understanding of 
their formation in the early universe. 

The investigation of the compact objects in the generalized scalar-tensor theories of gravity helps 
us understand them better. On the other hand, the investigation of matter in extreme conditions like 
these in the neutron stars may demonstrate new phenomena and new features of specific scalar-tensor 
theories of gravity, originating from the low energy limit of the string theory. Thus, at first time we 
may be able to reach theoretical indications of physical manifestation of the string theory in the real 
world H]. 

In the present paper we develop a direct numerical method for solving the equations of the general 
scalar-tensor theories of gravity including a dilaton potential term for the general case of mixed boson- 
fermion star. 

The physical motivation for considering mixed boson-fermion stars is connected with the fact that 
many of the present-day existing stars are of primordial origin being formed from an original gas of 
fermions and bosons in the early universe. That is why it should be expected that they are a mixture of 
both fermions and bosons in different proportions. The study of such mixed objects is a new interesting 
problem, whose investigation was started in reference There exist different candidates for boson 
fields in stars such as Higgs field of Standard model, or axion field being a pseudoscalar partner of 
dilaton in the superstring theory. They are unavoidable part of modern physics, nevertheless up to 
now we have no experimental evidence for their existence. Taking into account that according to the 
modern understanding of the initial state of universe a significant amount of these fields must have 
been present during the Big Bang phase, one has to expect some part of these fields to be present in 
the stars of primordial origin. The study of new observable effects of boson fields in such mixed stars 
may give new ways for discovery of the existence of the above hypothetical fields, which at present are 
the most intriguing new objects in modern physics. 

In the Einstein frame the field equations in the presence of fermion and boson matter are: 

Gl^nAn +Ti\+ 25.^9V -d'^di^Sl -f \u{^5l 
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V,VV+4t/'(¥') =-^o.^)(T + Tj 



(1) 



where Vi is the Levi-Civita connection with respect to the metric gij, [i — 0, ...,3; j = 0, ...,3). The 
constant is given by = SttG*, where G* is the bare Newtonian gravitational constant. The 
physical gravitational "constant" is G,A^((y9), where A{if) is a function of the dilaton field Lp depending 
on the concrete scalar-tensor theory of gravity. For example, in the framework of the Brans-Dicke model 
we have A(tf) = exp(-^7=^==), where ujbd is a parameter. 



The dilaton potential U{(p) can be written in the form U{(p) = m.'jjV{(p), where m^) is the dilaton 
mass and V{ip) is a dimcnsionlcss model function of (p. 

The complex scalar field 'I' describes the bosonic matter, while is its complex conjugated 
function. The quantity VF(*I'+\I/) is the potential of boson field, which can be chosen in the following 
form: ^ 

W{^+'^) = * - (*+^')^ 

where A is a parameter. 

The scalar function a{ip) = [In A{ip)\ determines the strength of the coupling between the dilaton 
field ip and matter. 

B F 

The quantities T and T arc correspondingly the trace of the energy-momentum tensor of the 
F B 

fermionic matter and the bosonic matter Tf . We note that in the present article we consider the 

fermionic matter only in macroscopic approximation, i.e., after averaging quantum fluctuations of the 

corresponding fcrmion fields. Thus, we actually consider standard classical relativistic matter. 

The explicit forms of the mentioned tensors are correspondingly; 

^- 1 

~A^{^) [a^^+a'* - 2A^{ip)W{^+-^)\^ 51 , (2) 

F 

T/= {e+p)uiu^ -p5i . (3) 

Here, the energy density and the pressure of the fermionic fluid in the Einstein frame are e = A'^{(fi)i 
and p = A'^{(fi)p, where e and p arc the physical energy density and pressure. Instead of giving the 
equation of state of the fermionic matter in the form p = p{e), it is more convenient to write it in a 
parametric form: 

e = eog{fi) p = eof{n), (4) 

where £o is a properly chosen dimensional constant, /x is the dimensionless Fermi momentum, and /(/x) 
and 5f(/x) are given functions (see below). 

The physical four- velocity of the fermionic fiuid is denoted by Ui. 

The field equations together with the Bianchi identities lead to the local conservation law of the 
energy-momentum of matter: 

^. F 

VjT^=a{^)Tdi^. (5) 

From now on, wc will take into consideration a static and spherically symmetric mixed boson- 
fermion star in asymptotic flat space-time. This means that the metric gij has the form: 

where r, 6, (f> arc usual spherical coordinates. 

The field configuration is static when the boson field 5* satisfies the condition: 

* = CT(r) e*"*. 

Here, u is a real number and a{r) is a real function. 

Taking into account the above-stated assumption, the system of the field equations is reduced to 
a system of ordinary differential equations (ODEs). Before writing the system explicitly, we are going 
to introduce a rescaled (dimensionless) radial coordinate by r ^ rnsr, r € [0, 00), where is the 
mass of the bosons (a prime will denote the differentiation with respect to the dimensionless radial 
coordinate r). 

We also define the following dimensionless quantities by: 



The components of the energy-momentum tensors of the fermionic and bosonic matter, written by the 
dimensionless quantities, are correspondingly: 



TS =bA\^)g{f,), Tl = Ti^ ^bA\^)f{ii), (7) 

Tl = \^^A\^) e-^a\T) + \a\^) V'^ - A\^)W{a\ (8) 

Tl ^-\^^A\^)e-^a\r)~-\A\^)e-^a'^-A\^)W{o''\ (9) 

Tl =-\^^A\^)e-^a\v)-r\A\^)e-^o'^ -A\^)W{o-'). (10) 



The parameter h — K^eQ/rn^ describes the relation between the Compton length of dilaton and the 
usual radius of neutron star in general relativity. 

It is necessary to note that two physically interesting borderline cases of pure bosonic and pure 
fermionic stars are formally contained in the above general system (0). For example, the model of pure 

F 

bosonic stars can be obtained from (|l|) by letting the tensor T/ to be zero. While the pure fermionic 
stars correspond to the field = 0. The case of pure bosonic stars in the scalar-tensor theories of 
gravity with a massive dilaton has already been discussed in our recent paper pst . In the present 
paper we consider the mixed boson-fermion stars. 



2 Formulation of the Problem 

Under the physical assumptions we have made, the field equations (|^) can be reduced to a system of 
ODEs. From mathematical point of view it is more convenient all ODEs to be of second order. That 
is why we first solve the Einstein equation G\ for e^: 



1 -r2 



Tl +ljW{ip) - i02A2((^)e-''cr2 - A'i{ip)W{a^) 



as a function of the quantities ^{x), v'{x), <y{x), <j'{x), (p{x), f'^x), and the spectral parameter f2, 
and then substitute the above expression in the other Einstein equations. In this way, in terms of the 
dimensionless quantities, the system of the field equations (Q) is reduced to the following system of 
ODEs: 



r 



+ — = 



+ (^0 - ^0 -2 T| + Tj? - Tl -2 T 



-l'V{v) + —bigl{T,+j'V{^)) 

ip' a{ip)/F 
-Z. + T + T 



1 



+ ^7V'M + ^(Ti+7V(^)) 
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r 



2A\ip)W'{a^)a 



{T,+j'Viip)) 



In the above equations, the potential of the bosonic matter W has the form: 

1 / . 1 



a 



-Aa 



and we suppose that W'{a^) 



Similarly, we set V'{ip) = 



2 V 2 



(11) 



(12) 



(13) 



The quantity Ti depends on the components of the energy-momentum tensors of the fermionic and 
bosonic matter (0)-(|9|): 

F F B B 

Ti = To" + T 1 + To" + Tl . 

B F 

The quantities T and T represent the traces of the these tensors, and are defined by the formulae: 

Correspondingly, the conservation law (|^) can be expressed as: 

.9(m) + /(m) 



The fermionic matter functions f{^) and g{p), entered in the above relations, have the form: 
fili) =i [(2Ai-3)VM + M2 + 31n (yTi + .JTni 
9{f^) =1 [(6Ai + 3)VM + M' -31n + ^TT, 



(14) 

(15) 
(16) 

Let us now complete the problem by adding proper boundary conditions (BCs) to the system of 
differential equations (11)-(^4|). 

The asymptotic flatness means that the function v{r) — > when r — s- oo. On the other hand, the 
nonsingularity condition at the center of the star requires the derivative i''(0) — 0. The same condition 
in relation to the dilaton field Lp{r) implies that the derivative <f'{0) = 0. At the same time, the 
function ip{r) at the asymptotic infinity (r — > oo) must be tpac = as it is required by the asymptotic 
flatness. The nonsingularity of the bosonic density (T(r) at the center of the star requires the derivative 
cr'(O) = 0. We need finite mass for the star, which implies (T(r) when r — > oo . In addition, the 
central value CTc = ct(0) must be given. Concerning the fermionic fiuid, we have to give the central 
density ic = s{0) or, equivalently, the central value /ic — /i(0). 

It should be noted that for the physically relevant equation of state of the fermionic matter there 
must be a point r = Rg < oo, where the pressure of the fermionic matter vanishes, i.e., Rg is the 
radius of the fermionic part of the star. 

As a conclusion, from the above-mentioned physical assumptions, we can formulate the following 
linear boundary conditions (BCs) for the quantities under consideration: 



^.'(0) = 0, 
^'(0) = 0, 
a'(0) = 0, 
fi{0) ^ He 



iy{(X)) ~ 0; 
</j(oo) = 0; 
cr(oo) = 0; 



(17) 
(18) 
(19) 
(20) 



Here, we denote (^(oo) hmr^oo(-)(^)- 

Apart from the unknown functions v(r\ (T(r), '~p(r\ and /i(r), the equation s ([Tl| )-(14) also include 
two unknown real parameters, i?s > and f2. However, the seven BCs (|l^-(|^are insufficient for 
their computation. In order to determine these parameters, we have to use additional conditions. 
In other words, the problem may be considered as a nonlinear eigenvalue problem, where R^ and 
are considered as "eigenvalues". For this purpose, further on we use two physically-clear additional 
conditions. 

The first one, given by the relation: 

(7(0) - a, (21) 

determines the density (Tc > of the bosonic matter in the star's center. The second one: 

ii{Rs) = 0, < i?^ < oo (22) 

describes the condition that the density of the fermionic matter must vanish at the radius of the star. 

Finally, we note that all the functions vir), (t{t), and (p{r) are defined in the whole real half-line 
r € [0,oo). It is easy to see that these functions are smooth in this interval including the point x ^ 1. 
Whereas, the fermionic density /i(r) is defined and smooth only inside the star, i.e., r € [0, i?s]. 



3 Method of Solution 



For solving the above formulated nonlinear eigenvalue problem the Continuous Analogue of Newton 
Method (CANM) (see - and comprehensive surveys |^) is applied. For convenience, a 
brief description of CANM can be found in the attached Appendix. 

The presence of the a priori unknown quantity Rs, however, is an obstacle for the direct use of 
CANM - the problem is the unknown internal boundary Rg. In order to overcome this obstacle, we 
introduce a new scaled coordinate x — r/Rg. As a result, the physical domain r e [0, oo) renders to 
the domain x € [0, cxd), and the star's radius r = Rg maps into the fixed point x = \. Then the BC 
(p^) for becomes 



m(i) = o. 



(23) 



Let Xi and X2 be two arbitrary points in the internal domain [0, 1]. We note that for the arbitrary 
functions /(/i), (/(/i), and a{Lp) the equation (|l4|) has a first integral, which can be presented as: 



/(m)+.9(a^) 



- {v2 — vi) + In 



where i^i, z/2, (pi, </?2i Mi: stand for the functions v{x)^ ^{x), fJ-{x) at the points xi and X2, respectively. 
Thus, for the model of the fcrmionic matter described by the conditions (p^, ( p^ we simply get the 
following algebraic equation: 



In 



(l+M2)A^(y2) 
(l+Ml)A2(^l) 



1^2 — i^l = 0. 



(24) 



For convenience, we introduce the vector y(a;) = {i^ix), ip{x), (t(x)}. Then the first three equations 
0) ~ (0) °^ ^^'^ problem and the corresponding BCs ( p7| ) - (|l9|) can be rewritten as follows: 



xy 



0, 



(25) 



y'(0) = 0, y((X3)-0, (26) 

where F — F{x, y, y', /i, Rg, fl) is 3D vector consisting of the right-hand sides (RHSs) of the equations 
(pl])-(p^) multiplied by R'^x. The differentiation with respect to the new independent variable x is 
denoted by (.)'. In the linear case, the advantages of such representation of the radial operator are 
discussed in p6[ . 

Following CANM, we introduce a "time-like" parameter t e [0, oo) and assume the unknown quan- 
tities depend on t as well: y = y{x, t), Rg ~ Rg{t), = Q{t). Let us suppose that the function /i = /i(a;) 
is known (see below). Then the CANM equations |Q corresponding to ( p5| ) and ( [2^ ) become: 

fdF \ , dF f 2 ^ dF\ dF 

= xy" + y'-F. (27) 
z'(0) = -y'(0), z((X)) = -y(^), (28) 

where E is an identity 3x3 matrix and 

y = z, Rg=p, n = to. (29) 



The respective Frechet derivatives at the point (y, Rg, fl) are dF/d{.) and the dot in ( |29| ) and below 
denotes the differentiation with respect to "time" t. 

The solution z{x) of the above equation is sought as a linear function towards the derivatives p 
and Lu 



z = u + pv + CJW, 



(30) 



where u{x),v{x), and w(x) are supposed to be new unknown 3D vector-functions of x. Substituting 
for them in equation (^), we obtain the following three vector ODEs of second order with respect to 
these quantities: 

dF dF 

-xu"-u' + |^u+g;u' =xy" + y-F (31) 
, dF OF , f 2 ^ dF\ 

,/ , 9F dF , dF 
-.w"-w' + -w + -w' (33) 

The above three equations are coupled with the following six BCs: 

u'(0) = -y'(0) , u(oo) = -y(^); (34) 
v'(0) = , v(oo) = 0; (35) 
vf'iO) = , vir(oo) = 0, (36) 



which are obtained from BCs (|2^), substituting for them with decomposition ( |30| ) also. Let us em- 
phasize that the above equations (pl|)-(p6|) have equivalent structures of the left-hand sides, which 
essentially facilitates their numerical treatment. 

In order to calculate the derivatives p and uj, we apply CANM for the first additional BC (pT|). 
This gives: 

a(0) - a, - (7(0). 

One more condition is required. Unfortunately, the second additional condition ( p3|) is not convenient 
for this purpose because knowledge about decomposition (|3^) concerning the function /i(x) is not 
available. We avoid this difficulty using the integral (|2j) for xi = and X2 = 1- Taking into account 
conditions ( [20| ) and (|23|), we obtain an algebraic equation with respect to the quantities i'(O), J^(l), 
ip{0), <p(l). After applying CANM to this equation, we get 

ln(l + Mc) - KO)] - 2 In ^Mjj = 0, 

where the abbreviation A' denotes the derivative of the function A with respect to the argument ip. 

Let us now eliminate all the derivatives in relation to "time" t by means of decomposition (^0|). 
As a result, we receive the following linear system of algebraic equations: 

aip + biLO = ci 

(37) 

with respect to the unknown derivatives p and uj. The coefficients in formulae (^7|) are given by: 

c, ^ m [1 + „j - Ml) - Ko)l - 2^ u,m + „,,„) 

02 = ^3(0), = 1^^3(0), C2 = 0-c - cr(0) - U3(0). 



Obviously, the explicit form of the coefficients in system (|37| ) depends on the concrete choice of 
functions f{p) and g{p). 



4 General Sequence of the Algorithm 



We discretize the continuous "time-like" parameter t G [0,oo) in the following way: t^+i = tk + t^., 
to = 0, where k = 0,1,2,... denotes the number of iterations, and the "time" step Tk is generally 
assumed as a variable quantity. Next, we use the Euler difference scheme [pll to approximate the 



"time" derivatives in equations (29). Then we can write 



yfe+i (x) = yfe (x) + Tk [ufe (x) + pkVk (x) + ujkWk (x)] , 
Rs,k+1 =Rs,k + TkPk, (38) 

Let us suppose that the functions Vkix), fkix), (Jk{x), pk{x) and the parameters i?s,fe, ^k are given. 
We solve the linear BVP (|3l])-(p3[) and, thus, we compute the functions Ufc(a;), Vfe(a;), •Wk{x). Next, 
to obtain the derivatives pk and ojk we solve system (^. After that, using decomposition ( |38| ) for a 
selected tj,, we calculate the functions Vk+i(x), ipk+i{x)^ Uk+i(x), the radius of the star Rg^k+i, and 
the quantity flk+i as well at the new stage k + 1. In the end, we calculate the function p,k+i{x) at 
the new stage, according to the recurrent formula, which can be obtained immediately from the first 
integral (H). 

For every iteration k an optimal time step Topt is determined in accordance to the Kalitkin& 
Ermakov formula , : 

~ 6(0) + s{iy 

where the residual S{t) is represented as follows: 

5{Tk) = max [6f, {Rs,k + TkPkf , i^k + TkUJkf] 



and Sf is the Euclidean residual of RHS of the equation (|3l|). Formula (|39[) provides approximately 
the minimal value of the residual for the current solution, given by (^8|). 

The criterion for termination of the iterations is S{Topt) < £, where e 10^^ 10^^^. Then, for 
the sought solutions we set i^{x) = Vk+iix), ip{x) = ipk+i{x), a{x) = ak+i{x), = i?s,fe+i, ^ = ^k+i- 

The use of the standard programs available, for example, via the Internet to solve numerically 
the linear BVPs (|3l|)-(|3^) is unhandy for many reasons. Because of that, the spline-collocation scheme 
is employed in our case. 

We introduce a nonuniform grid 

A : Xi+i = + /li, i = 0,1, . . . , Ns, Ns+i, . . . , N - 1, xq = 0, xn = X^o, 

on the interval x G [0,Xoo], condensing to the points x = Q and x = \. Here, Xoa is the "actual 
infinity", Ns is the number of the node x = 1, is the full number of the subintervals, and hi is the 
grid step. We will seek approximate solutions of the above linear BVPs as a cubic spline on the grid 
A. Namely, for x £ [xi, Xi+i], i = 0, N ~ 1 we set 

V{x) = MO) U. + MO) M, + MO) U,;+i + MO) M,+i. (40) 

In the above formula the relative coordinate — {x-'Xi)/hi and the known functions ipi{0), I — 1, . . . , 4, 
are the coefficients of the spline. For simplicity in the last formula, we introduced the 3x3 matrices 
U and M, consisting of the coordinates of the vectors u, v, w from ( ^ ) and their first moments at 
the spline nodes Xi , i — 0, . . . , N . According to the collocation method in every subinterval 
[xi,Xi-^-i] , i = 0, A — 1, the system ( |3T| ) - (33) is satisfied at the corresponding Gaussian points 



9i = 1/2 — and 62 = 1/2-1- •\/3/6. This kind of discretization yields an algebraic system with 

respect to the functions and their moments at the spline nodes. The corresponding matrix has an almost 
block-diagonal structure (see [Q). Therefore, at the i— th block (i ~ 1, . . . , A^ — 1) the collocation 
equations have the form: 



4J1 \\blJ HJ KJ 



\bU llcLII 114. 



U,+i 
V M,+i J 



where is the vector of RHSs of the equations (|31|)-(|33|) at the coUocation nodes, while the superscript 
corresponds to the number of these nodes. Here: 



^ " + + (^)., ^^^^ + 

for fc= 1,2,3, 71=1,2,3, j = l,2, 

and the quantities ^ij — Xi + 9jhi are the absolute coordinates of the collocation points. The derivatives 
of the spline coefficients ipi with respect to the relative coordinate 9 are dotted. 

The dimensions of the first and the last blocks in the global matrix are greater, since we add two 
matrix rows corresponding, respectively, to the left and right BCs. 

Formula (^0|) is also used for the approximation of the RHSs of system (^l]) - (|3^ ) in the collocation 
points. 

The spline-difference schemes of this kind have a high order of approximation 0(h'^) , where h = 
max{ft.i}, i = 0, N. 

It is clear that for solving all the three algebraic systems, corresponding to the linear BVPs (|l|) - 

at every iteration, only one L [/-decomposition is necessary. 
Depending on the initial values of the government physical parameters, the number of iterations 
varies approximately in the range 4 ~ 16. If we vary some solution as a function of one of the parameters 
Mci fcj 7j a, or 6, then we use the previous solution as an initial approximation for computing the next 
one. 



5 Results and Discussion 

In order to be specific in the present article, we focus our attention on a concrete scalar-tensor gravity 
model, characterized by the functions 

AM=exp(^) and V{^) ^ (1 - [A{ip)]-'f . 

For more details concerning this gravitational model, we refer the reader to the recent paper ||39[| and 
the references therein. 

The order of approximation of the used spline-difference scheme is verified by the Runge rule. 
The Runge rule is presented by the formula: 

where p is Runge's number and yu, yh, yn are the values of the grid function y at the given node, 

2 4 

computed on meshes with steps h, ft./2, and h/A. In our case p must be approximately equal to 4. 

In Table 1 the values of the sought grid functions at the point x — 1, the corresponding radius of 
star and the quantity Vl for CTc = 0.8, /ic = 1, A = 0.01, 7=1,6=1, and Xoo — 128, are shown. 

Therefore, it is obvious that the Runge relationship is satisfied both for the functions and the 
eigenvalues Rs and Vl. 

The correctness of the spline-difference scheme is verified through appropriate numerical experi- 
ments consisting of both grid doubling and doubling of the "actual infinity" . For this purpose, uniform 
meshes are used with numbers of the spline nodes N = 256, 512, 1024, 2048, respectively. It turns 
out that the relative error between the values of the functions v{x), (p{x), and (t{x), varies in the 
range 0.1% - 1% when the mesh is "coarse" {N = 256, 512), and in the range 0.003% - 0.02% when 
the mesh is "fine" {N = 1024, 2048). Similar experiments are carried out with the "actual infinity" 



Table 1: Data for checking the Runge rule. 



h 








Rs 


n 


i 


-1.0059230404 


-0.0471137759 


0.4777335163 


1.1609111685 


0.8006662485 


64 


-1.0059334054 


-0.0471120738 


0.4777483180 


1.1608888836 


0.8006671950 


-1.0059342032 


-0.0471119781 


0.4777490917 


1.1608875328 


0.8006672467 


P 


3.61 


4.22 


4.37 


4.06 


4.28 



Xoo = 64, 128, 256. It is interesting to note that the relative error between the set functions if{x) and 
a{x) is very small (less than 10~'*%), while the function i'{x) is more sensitive with respect to the 
choice of the quantity Xoo- This fact is fully explainable if we take into account that the function v^x) 
decreases slowly at the infinity compared to the other functions. (Theoretically z^(x) ^ when 
a; — > oo. Here, the quantity M is the total star mass.) The computed values of the derivative h''{Xao) 
as a function of the "actual infinity" Xoc are presented in Table 2. It is easy to see the relationship 
i>'{Xao) — where the constant C > depends on the concrete solution (for the above solution 
C w 1.133). 



Table 2: Asymptotic behaviour of the derivative i^' at the "actual infinity" X, 





32 


64 


128 


256 


512 


ly'iXoo) 


1.07246x10"^ 


2.63721x10"'' 


6.53945x10"'^ 


1.62825x10"'^ 


4.06241x10"" 



All government parameters are varied in wide physically-admissible ranges. As initial distributions 
of the functions v{x), ip{x), cr{x) and fJ.{x) both analytic and numerical approximations are used. 

Results concerning a family of solutions will be considered below. They are obtained for the 
following fixed values of the parameters: /ic — 0.5, A = 10, 7 = 10, 6 = 1, and the "actual infinity" 
Xoo — 128, when the parameter CTc runs the interval [0.1,0.9]. 
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Figure 1: The function v{x) in dependence on the parameter ac'- "o" - ctc = 0.1; "A" - ac — 0.5; "V" 
- Ctc 0.9. 

Figure 1 presents the dependence of the function i'{x) on the dimensionless coordinate x for three 
different values of the central bosonic density CTc. It is seen that when CTc increases, the absolute value 
of i'{x) as a whole decreases and at great distances (from 3 star radii when CTc = 0.1 until 45 star radii 
in the case CTc = 0.9) from the star's center approaches asymptotically zero. The qualitative behaviour 
of the three curves, however, remains the same. Such a behaviour is natural and should be expected if 
the differential equation ( |ll|) for v(r) is taken into account. From physical point of view, this behaviour 
is natural also because the function exp(^^^) is related to the gravitational potential. 




Figure 2: The potential of dilaton (p{x) as function of the parameter ac- "o" - ctc = 0.1; "A" - ac — 0.5; 
"□" - ac = 0.7; "V" - CTc = 0.9. 

Figure 2 presents the dependence of the dilaton field f{x) on the dimensionless coordinate x for four 
different values of CTc- The qualitative behaviour of the field f{x) as a function of a^. is the following. 
For small values when a^. increases, the dilaton field around the center of the star decreases. Then, 
after some critical value a* the behaviour of ip{x) is changed and ip{x) around the center of the star 
begins to increase with the increase of (Tc. The cause of the described behaviour is the presence of 

B J—. B 

the term T on the RHS of equation (p_2p. For sufhciently small values of the density ac the term T is 

negative and has a dominant contribution with respect to the term T- For the sufficiently large central 

B ^ ^ ^ . . P . 

value (Tc (cTc > f *), the term T changes its sign and amplifies the contribution of T, leading to the 

increase of the function ^{x). 



Parameters: = 0.5, A = 10, y = 10, = 1 
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Figure 3: The bosonic density a{x) in dependence on the parameter ac- "o" - tXc = 0.1; "A" - ac — 0.5; 
"V" - ac = 0.9. 

From a physical point of view, the described behaviour of the dilaton field (and consequently the 
behaviour of the physical gravitational "constant" G'*A^((p)) for the central values ac > a* seems to be 
strange. In order to clarify this situation, we have to take into account that in the range ac > a* (for the 



fixed value of the central fermionic density fic) the star is unstable and, therefore, the mentioned range 
is not physically relevant. Such a behaviour has to be considered only as an iteresting mathematical 
fact. In the domain of stability < CTc < cr*, as we have already seen, the dilaton field ip{x) has a 
normal physical behaviour - it decreases when the parameter ctc increases. 



Parameters: = 0.5, A = 10, Y= 10, ^ = 1 




Dimensionless coordinate x 



Figure 4: The fermionic density /i(a;) in dependence on the parameter ac- "°" - '^c = 0.1; "A" - 
fJc 0.5; "V" - (Tc = 0.9. 

The dependence of the bosonic density cr{x) on the dimensionless coordinate x for three different 
values of <Tc is presented on Figure 3. The qualitative behaviour is the same for all three different 
values of CTc- It approaches zero at infinity (rapidly when CTc = 0.1 and more slowly when Cc increases). 
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Figure 5: The radius of the star and the quantity f2exp( ) as functions of the parameter ctc- 




In the next Figure 4 the dependence of the fermionic density fj,{x) on the dimensionless coordinate 
X is presented for three different values of CTc- The qualitative behaviour of the three curves is similar. 
In agreement with the initial assumption, it is nontrivial only within the star. It is seen that when the 
value of (Jc increases, the density ^{x) increases as a whole, too. This fact is related to the effect of an 
increase of the gravitational field with the increase of (Tc - the star becomes more compact, which leads 



to the greater density of matter, respectively to the function ^i{x). The same may be seen in Figure 5 
- when the central value CTc increases, the radius of the star Rg decreases about 10 times. 

From physical point of view it is important to get knowledge about the behaviour of the quantity 
as a function of the central value tic- That quantity may be considered as the energy of 
one boson particle in the gravitational field yielded by the rest matter (in the Einstein frame) . Figure 
5 clearly shows that the quantity f2exp(— ■^^^) increases along with Uc. Such behaviour should be 
expected, because the energy of the system has to increase with the increase of the central density CTc 
of the star. 

Concluding Remarks 

Based on CANM an iteration method for solving the nonlinear BVP, describing a static spherically- 
symmetric boson-fermion star, is developed. 

A linearization of the main equations of the star renders the original two-parametric nonlinear 
spectral problem to three two-point linear vector BVPs and a linear system of algebraic equations 
for the spectral parameters (the radius of the star Rg and the frequency of the bosonic field). A 
spline-collocation scheme of fourth order of approximation for solving numerically these BVPs is used. 

Our basic physical result is that the structure and the properties of the star in the presence of a 
massive dilaton field depend essentially on both its fermionic and bosonic components. This shows 
that a careful investigation of these properties may give new real ways to discover physical effects of 
the hypothetical boson fields and dilaton field in stars. 

Appendix 

For reader's convenience, we briefly explain the main ideas of CANM. 

CANM can be treated as a particular case of the continuous analogues of iteration methods, strictly 
formulated and studied by M.K. Gavurin in 1958 (see the review in |4^). Among the number of papers 
devoted to the theoretical development and applications of CANM for solving wide classes of nonlinear 
equations, we will indicate the basic papers [|9| - |Q as well as the reviews @, 

Let us consider the nonlinear equation: 

x{y) = 0, (41) 

where xiv) is an operator defined in a Banach space Y. We suppose that the equation ( ^l| ) has an 
isolated exact solution y* g Y. Let the element j/o G Y (an initial approximation to y*) be given. To 
solve equation (pi]), we can use an iteration process, usually taking it in the form: 

yn+i^yn + i^ivn), n = 0,1,2,.... 

Here, n indicates the number of iterations and ip is an appropriate function, which carries Y into itself 
and has the same zeroes as x- 

The choice of the function ip{y) depends on the kind of concrete iteration method used. 

According to Gavurin's idea, for each iteration process of such kind it is possible to formulate the 
corresponding continuous analogue in the following way. Let us consider an abstract function y(t) of 
the independent continuous variable t e [0,oo) instead of the sequence yo,yi, ■■■,yn, and suppose 
that y{tn) = yn for each n. Then, we can introduce the derivative y{t) instead of the increment 
yn+i — yn and replace (^) with the abstract initial value problem on the interval t G [0,oo) 

m^Hv), y(0) = yo. (42) 

Such a transition from a difference equation to a differential one has many advantages, both in 
pure theoretical and applied aspects. 

In the case of Newton's method, we set ip{y) — — x'(!/)^^x(y)j where x'iy) is the corresponding 
Frechet derivative of xiv)- Then, the main equation of CANM, arising from (^), can be rewritten in 
the form: 



x'{y)y = ~x{y)- 



(43) 



Obviously, the above ODE has a significant first integral of the kind: 



x(2/W) = x(yo)e-*, (44) 

which means that xiui^)) ~^ when < — > cxd. 

Various theorems, based on concerning the convergence of a path y{t) to the exact solution 
y* have been proved. For example, a theorem jsj], which guarantees the convergence of CANM for a 
simple BVP, is cited below. 

The following BVP is considered: 

-y" + f{x,y)=0, xe(0,l), (45) 
y(0) = 0, y(l) = 0. (46) 



Theorem 1 Let the BVP ^4^, have an isolated solution y*{x) and: 

i) the function /(x, y) have continuous partial derivatives up to the second order in some domain D; 

ii) the linear BVP 

-w" + f'y{x,y)w^Q, xe(0,l), 
w(0) = 0, w(l) = 0, 

have only a trivial solution w{x) = for every smooth function y{x) € D; 

iii) the initial approximation yo (x) £ D he a smooth enough function satisfying: 

\\-Vq+ f{x,ya\\<e for e>0. 



Then the system 



W + f'{x,y)w = y" - f{x,y), y = w, 



with BCs w{0,t) — 0, w{l,t) — 0, and an initial condition y(x,0) = yoix), has in [0,1] U [0, oo) an 
unique solution, satisfying the relation: 

lim \\y{x,t) - y*{x)\\c2[Q,i] =0. 

The numerical solution of CANM equation ( ^ ) is based on an appropriate scheme for discretization, 
which has to be stable for the asymptotic stability of the path y{t). The most frequently used one is 
Euler's scheme (see the details in the above cited papers). At first, the linearized equation: 

X'{yn)Wn = -X(yri), (47) 

is solved with respect to the increment Wn, and then the next approximation is obtained via the 
formula: 

2/«+l =yn+ TnWn- (48) 

Here, < t„ < 1 is an iteration parameter. When t„ = 1, the classical Newton method is obtained. 
We note that the choice of t„ is important for the rapid convergence of the process. It is possible to 
choose this parameter so that the range of convergence is wider in comparison to the classical Newton's 
method ||. 

Theorems regarding the convergence of iterations (^7|), ( p8[ ) for wide enough hypotheses as well as 
essential generalizations of CANM, are discussed in the above cited papers. 
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